clear all
*cap log close
set more off
set seed 603

// other stuff -------------
local nbin = 16
local Q    = 2
// -----------------------------------------------------


// --------------------------------------------------------------------------------
// quadratic extrapolation figure 

* read in polynomial estimates ------
import delimited using "${est}/recid_poly.csv", clear 
keep if poly == `Q'
replace par = trim(upper(par))

* store estimates for printing -------
qui summ est if par == "Y0"
local b = r(mean)
qui summ se  if par == "Y0"
local s = r(mean)
local py0 = "E(Y0) = `:di %4.3f `b'' (`:di %4.3f `s'')"

qui summ est if par == "Y1"
local b = r(mean)
qui summ se  if par == "Y1"
local s = r(mean)
local py1 = "E(Y1) = `:di %4.3f `b'' (`:di %4.3f `s'')"

* store estimates for graphing ---------
keep if par == "Y0"|par == "Y1"
gen z = (par=="Y1")
gen flag = 1 
gen y    = est
keep z y ub lb flag 
tempfile est 
save    `est'

* store binscatter ------------------
use "${data}/out/4-main", clear 
gen Y = cite_ny1 
keep Y Z totfe 
binsreg Y Z, absorb(totfe) binspos(es) nbins(`nbin') savedata(bs) replace 

* store polynomial fit ----------------
forval q = 1/`Q' {
	gen z_`q' = Z^`q'
}
reghdfe Y z_*, absorb(totfe) 
local b0 = _b[_cons]
forval q = 1/`Q' {
	local b`q' = _b[z_`q']
}

* setup data for graphing ------------
use bs, clear 
rename dots_x   z 
rename dots_fit y
qui append using `est'
gen qfit = `b0'
forval q = 1/`Q' {
	replace qfit = qfit + `b`q''*z^`q'
}
rm bs.dta 

* create graph -----------------------
sort z
#delimit ;
scatter y z, msymbol(Oh) mcolor(dknavy) msize(medium) || 
line qfit z, lcolor(dkgreen) ||
scatter y z if flag==1, mcolor(dkgreen) msize(med) ||
rspike ub lb z if flag==1, lcolor(dkgreen)
graphregion(color(white)) plotregion(lcolor(black) lwidth(medthin)) 
ylab(0.33(0.01)0.38,nogrid) xlab(,nogrid)
xtitle("Officer Stringency") ytitle("Pr(Reoffend)")
legend(ring(0) cols(1) pos(2) order(1 2) region(lstyle(border))
lab(1 "Local Mean") lab(2 "Quadratic Fit"))
text(0.332 0.2 "Extrapolation Estimates:" "`py0'" "`py1'", place(n) size(medium)) ;
#delimit cr
graph export "${out}/main/baseline_extrap.pdf", replace 
// --------------------------------------------------------------------------------



// --------------------------------------------------------------------------------
// levels and gains plot 

* read in polynomial estimates ------------
import delimited using "${est}/recid_poly.csv", clear 
keep if poly == `Q'
replace par = trim(upper(par))

* store info for printing ----------------
qui summ est if par == "Y0D1" 
local b = r(mean)
qui summ se  if par == "Y0D1"
local s = r(mean)
local pt1 = "E(Y0) = `:di %4.3f `b'' (`:di %4.3f `s'')"

qui summ est if par == "ATT" 
local b = r(mean)
qui summ se  if par == "ATT"
local s = r(mean)
local pt2 = "ATT = `:di %4.3f `b'' (`:di %4.3f `s'')"

qui summ est if par == "Y0D0" 
local b = r(mean)
qui summ se  if par == "Y0D0"
local s = r(mean)
local pu1 = "E(Y0) = `:di %4.3f `b'' (`:di %4.3f `s'')"

qui summ est if par == "ATU" 
local b = r(mean)
qui summ se  if par == "ATU"
local s = r(mean)
local pu2 = "ATU = `:di %4.3f `b'' (`:di %4.3f `s'')"

* populate levels and gains ----------------
do "${code}/utility/create_levels_gains.do"

drop if mi(y0)
replace ind = 0.2 if _n == 1
replace ind = 0.5  if _n == 2 
replace ind = 0.8 if _n == 3



* build plot -----------------------------
#delimit ; 
twoway 
(rbar y1u y1l ind, color(maroon%20) barwidth(0.02))
(scatter y0 ind, msize(large) msymbol(Oh) mcolor(navy))
(pcarrow y0 ind y1 ind, lcolor(maroon) lwidth(medthick) mcolor(maroon) barbsize(medlarge)),
xscale(r(0 1))
xlabel(0.2 "All Motorists" 0.5 "Treated" 0.8 "Untreated", nogrid) ylab(,nogrid) 
legend(ring(0) pos(2) order(2 3) cols(1) region(lstyle(border))
label(2 "Untreated Level") label(3 "Treatment Effect")) 
xtitle("") ytitle("Pr(Reoffend)") 
graphregion(color(white)) bgcolor(white) plotregion(lcolor(black) lwidth(medthin)) 
text(0.31 0.2 "Treated:" "`pt1'" "`pt2'", place(n) size(medium))
text(0.29 0.2 "Unreated:" "`pu1'" "`pu2'", place(n) size(medium));
#delimit cr
graph export "${out}/main/baseline_est.pdf", replace 
// --------------------------------------------------------------------------------



// --------------------------------------------------------------------------------
// robustness table ----------

* initialize LaTeX table ----------
capture erase "${out}/main/table_robustness.tex"
file open fh using "${out}/main/table_robustness.tex", write replace

* table header --------------------
file write fh ///	
	"&\multicolumn{4}{c}{\$Y_0\$} && \multicolumn{4}{c}{\$Y_1 - Y_0\$}\\" _n /// 
	"\cmidrule(l){2-5} \cmidrule(l){7-10}" _n ///
	"& (1) & (2) & (3) & (4) && (5) & (6) & (7) & (8)\\" _n ///
	"& All & \$D=1\$ & \$D=0\$ & Diff && All & \$D=1\$ & \$D=0\$ & Diff \\" _n ///
	"&     &         &         &      && (ATE) & (ATT) & (ATU) &  \\" _n ///
	"\hline" _n

* first panel: polynomials -----------
file write fh ///
	"\multicolumn{6}{l}{\textit{\underline{Panel A: Polynomials}}}  \\"  _n ///

import delimited using "${est}/recid_poly.csv", clear
replace par = trim(upper(par))
local glab2 = "$ q=2$"
local glab3 = "$ q=3$"
local glab4 = "$ q=4$"
local glab8 = "$ q=8$"

foreach q in 2 3 4 8 {

	local label = "`glab`q''"
	local k = 1
	foreach p in Y0 Y0D1 Y0D0 Y0D1MY0D0 ATE ATT ATU ATTMATU {
		qui summ est if poly == `q' & par == "`p'" 
		local m`k' = "`:di %4.3f `r(mean)''"
		qui summ se  if poly == `q' & par == "`p'"
		local s`k' = "(`:di %4.3f `r(mean)'')"
		local ++k 
	}	

	*di "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\"
	*di "        & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 

	file write fh "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\" _n
	file write fh "          & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 
}	

* second panel: local means -------------
file write fh "[1em]" _n 
file write fh ///
	"\multicolumn{6}{l}{\textit{\underline{Panel B: Local means}}}  \\"  _n ///

import delimited using "${est}/recid_bins.csv", clear
keep if type == "mean"
replace par = trim(upper(par))
replace bw  = int(bw*100)
replace bw  = 1 if bw == 0

local glab1 = "$ bw=0.1$"
local glab2 = "$ bw=0.05$"
local glab3 = "$ bw=0.01$"

local j = 1
foreach q in 10 5 1 {

	local label = "`glab`j''"
	local k = 1
	foreach p in Y0 Y0D1 Y0D0 Y0D1MY0D0 ATE ATT ATU ATTMATU {
		qui summ est if round(bw,0.01) == `q' & par == "`p'" 
		local m`k' = "`:di %4.3f `r(mean)''"
		qui summ se  if round(bw,0.01) == `q' & par == "`p'"
		local s`k' = "(`:di %4.3f `r(mean)'')"
		local ++k 
	}

	*di "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\"
	*di "        & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 

	file write fh "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\" _n
	file write fh "          & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 
	local ++j
}	

* third panel: propensity score --------------
file write fh "[1em]" _n 
file write fh ///
	"\multicolumn{6}{l}{\textit{\underline{Panel C: Polyomials in propensity score}}}  \\"  _n ///
	
import delimited using "${est}/recid_pscore.csv", clear
replace par = trim(upper(par))
local glab2 = "$ q=2$"
local glab3 = "$ q=3$"
local glab4 = "$ q=4$"
local glab8 = "$ q=8$"

foreach q in 2 {

	local label = "`glab`q''"
	local k = 1
	foreach p in Y0 Y0D1 Y0D0 Y0D1MY0D0 ATE ATT ATU ATTMATU {
		qui summ est if poly == `q' & par == "`p'" 
		local m`k' = "`:di %4.3f `r(mean)''"
		qui summ se  if poly == `q' & par == "`p'"
		local s`k' = "(`:di %4.3f `r(mean)'')"
		local ++k 
	}

	*di "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\"
	*di "        & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 

	file write fh "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\" _n
	file write fh "          & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 
}	

* fourth panel: within locations --------------
file write fh "[1em]" _n 
file write fh ///
	"\multicolumn{6}{l}{\textit{\underline{Panel D: Witin-locations}}}  \\"  _n ///

local file1 = "recid_troop"
local file2 = "recid_county"
local file3 = "recid_spc"
forval j = 1/3 {
	import delimited using "${est}/`file`j''.csv", clear
	gen spec = `j'
	if `j' == 1 {
		tempfile temp
		save    `temp'
	}
	if `j' >= 2 {
		qui append using `temp'
		tempfile temp
		save    `temp'
	}
}

use `temp', clear 
replace par = trim(upper(par))
keep if poly == 2 
local glab1 = "Troops"
local glab2 = "Counties"
local glab3 = "FM (2022)"

foreach q in 1 2 3 {

	local label = "`glab`q''"
	local k = 1
	foreach p in Y0 Y0D1 Y0D0 Y0D1MY0D0 ATE ATT ATU ATTMATU {
		qui summ est if spec == `q' & par == "`p'" 
		local m`k' = "`:di %4.3f `r(mean)''"
		qui summ se  if spec == `q' & par == "`p'"
		local s`k' = "(`:di %4.3f `r(mean)'')"
		local ++k 
	}

	*di "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\"
	*di "        & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 

	file write fh "`label' & `m1' & `m2' & `m3' & `m4' && `m5' & `m6' & `m7' & `m8' \\" _n
	file write fh "          & `s1' & `s2' & `s3' & `s4' && `s5' & `s6' & `s7' & `s8' \\" _n 

}	


file write fh "\bottomrule" _n  

file close fh
macro drop fh
// --------------------------------------------------------------------------------






